# this should make a lizard walk around in some 3d matrix, collecting temperature
# the lizard walks randomly
# there is no regard for body temperature
# there is a fixed Pr(moving) = (0.5/12). this means the lizards has a 50% chance
# of moving after 2 minutes
# time steps are every 10s
# Tb is evaluated every time step

require(deSolve)
library(deSolve)



# *************************************************************************                               
#build environmental temperature array
Te_array <-function (Te)
{
  a<-length(Te)/100
  Temp<-array(Te, dim=c(a,10,10))
  return(Temp)
}
  
  
  

# **************************************************************************
# to find Tb. This needs library(deSolve)
calc_Tb <- function(Time, State, Pars)
{
  with(as.list(c(State, Pars)),
  {
    DT <- -k*(Tb_init-(Ti+(Tf-Ti)*t))
    return (list(DT))
  })
}                                  

 


# **************************************************************************       
# this figures out the probability of moving
           
Prmove <- function (Pr, Tb, VTMax, VTMin, TselUpper, TselLower)
{
  if (Tb > VTMax)
  {
    probmove<-1
    return(probmove)
  }
  if (Tb < VTMin)
  {
    probmove<-1
    return(probmove)
  }
  if (Tb > TselUpper)
  {
    #this calculates Prmove with probability from a straight line interpolation
    #starting at lowest (Tsel, pr/move) and ending at highest (VTMax,1)
    probmove <- rbinom(1,1,(((1-(Pr/move))/(VTMax-TselUpper))*(Tb-VTMax)+1))
    return(probmove)
  }
  if (Tb < TselLower)
  {
    #this calculates Prmove with probability from a straight line interpolation
    #starting at lowest (Tsel, pr/move) and ending at highest (VTMin,1)
    probmove<- rbinom(1,1,((((Pr/move)-1)/(TselLower-VTMin))*(Tb-VTMin)+1))
    return(probmove)
  }
  else probmove <- (rbinom(1,1,(Pr/move))) 
  return(probmove)
} 
              
              
              
              
# ************************************************************************** 
# this fingures out where the lizard steps in the moore neighbourhood 
    
steps <- function (probmove)     
{ 
  if (probmove==1)
  {
    #randomly pick a number in the moore-neighbourhood
    moore.step<-round(runif(1,min=1,max=8))
  }
  else
  {
    #new location = old location
    moore.step<-9
    #z=rbind(z,z[i,])
  } 
  return(moore.step)         
}
 
 



# **************************************************************************
#this one makes a torus so the lizard doesn't walk out of bounds
torus <- function (x, y)
{
   #I need a lookup table to turn moore-cell coordinates into cartesian coordinates
  # when standing at 9, the moore-neighbourhood looks like:
  #   1 2 3
  #   4 9 5 
  #   6 7 8
  
  lookup.cell<-array(c(-1,0,1,-1,1,-1,0,1,0,1,1,1,0,0,-1,-1,-1,0),dim=c(9,2))
  coord <- lookup.cell[x,]
  loc_coord <- y + c(0,coord)
  #plus1 <- array((loc_coord+c(0,1,1)), dim=c(1,3))
  #mod <- c(plus1[,1], plus1[,2]%%10, plus1[,3]%%10)
  #new_loc <- mod + c(0,-1,-1)
  new_loc <- array(c(loc_coord[1], loc_coord[2]%%10, loc_coord[3]%%10), dim = c(1,3))
  if (new_loc[,2] == 0)
  {
    new_loc[,2] = 10
  }
  if (new_loc[,3] == 0)
  {
    new_loc[,3] = 10
  }
  return(new_loc)
}
   
  
  
  


# **************************************************************************
# this is the function that takes a step, gets Te, works out Tb, find Pr(move) and takes another step

Lizard_move <- function (Te=Te, TselLower=28, TselUpper=32, VTMax=39, VTMin=15, CTMax=44, CTMin=4, Pr=0.5, move=12)
{     
  # find the point where the lizard starts, to give it a starting env_temp and Tb.
  start.cell<-array(c(round(runif(2,min=1,max=10))),dim=c(1,2))
  location<-array(c(1,start.cell[,1],start.cell[,2]),dim=c(1,3,1))
  Temp<-Te_array(Te)
  Tenv<-c(Temp[location[1],location[2],location[3]])
  Tb<-Tenv
  #Tb starts at environmental temperature
  

  while (dim(location)[1]<dim(Temp)[1])
  {    
    for (i in 1:60)
    {
      Prob.move <- Prmove(Tb[i], VTMax, VTMin, TselUpper, TselLower)
      step.move <- steps(Prob.move)
      #location<-rbind(location, tail(location, n=1)+(c(0,c(lookup.cell[step.move,]))))
      location <- rbind(location, torus(step.move,c(location[length(location)/3],location[length(location)*2/3],location[length(location)])))      
      Tenv <- c(Tenv, Temp[location[length(location)/3],location[length(location)*2/3],location[length(location)]])
      #Pars <- c(k = 0.6, Ti = 28, Tf = 29, t=1)
      Pars <- c(k = 0.6, Ti = Tenv[length(Tenv)-1], Tf = Tenv[length(Tenv)], t=1)
      yini <- c(Tb_init = Tb[length(Tb)])
      times <- seq(0, 10, by = 9)
      #Tbnew <- calc_Tb(Tenv[i],Tb[i])
      #Tbtemp <- ode(func = calc_Tb(Pars=Pars), y=yini, parms = Pars, times = times)
      Tbtemp <- ode(func = calc_Tb, y=yini, parms = Pars, times = times)
      Tb<-c(Tb,unname(Tbtemp[nrow(Tbtemp),ncol(Tbtemp)]))
    }
    #want to keep (location, Tenv, Tb)
    location<-rbind(location, tail(location,n=1)+(c(1,0,0)))
  }
  
      
#end of lizard_move           
}      
     
     
# to run it: 
#library(deSolve)
#source("lizard-moves-random-4.r")
#Lizard_move(Te, TselLower, TselUpper, VTMax, VTMin, CTMax, CTMin, Pr, move)              